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' Abstract. This talk reports on the status of an approach to the numerical study of 

, isolated systems with the conformal field equations. We first describe the algorithms 

used in a code which has been developed at the AEI in the last years, and discuss a 
?H ' milestone result obtained by Hubner. Then we present more recent results as examples 

to sketch the problems we face in the conformal approach to numerical relativity and 
outline a possible roadmap toward making this approach a practical tool. 
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1 Introduction 



Since the early work of Hahn and Lindquist |j| , SmarrQ and Eppley ||, numer- 
ical relativity has become an important subfield of the theory of gravitation. To 
outsiders the progress often seems marginal and unsatisfactory. The classic goal 
of providing waveform catalogs for the newly built gravitational wave detectors 
' has still not been reached (although considerable progress has been made recently 

, [|J). By the nature of general relativity, the simulation of isolated systems poses 

particularly hard problems. Mathematically such systems can be formalized by 
O the concept of asymptotically flat spacetimes (see e.g. the standard textbook 

of Wald ||), but it turns out that important quantities such as the total mass, 
?-h . (angular) momentum or emitted gravitational radiation can only consistently be 

defined at infinity. The traditional approach of introducing an arbitrary spatial 
cutoff introduces ambiguities and is not satisfactory at least from a mathematical 
point of view. A remedy is suggested by conformal compactification methods, 
such as the characteristic approach presented by Luis Lehner in this volume, 
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or Friedrich's conformal field equations, which he describes in this volume. The 
latter approach avoids the problems associated with the appearance of caustics 
in the characteristic formulation by allowing to foliate the compactified metric 
by spacelike hypersurfaces. These hypersurfaces are analogous to the standard 
hyperboloid in Minkowski spacetime and are asymptotically null in the physical 
spacetime. The price to pay is the loss of the simplicity inherent in the use of 
null coordinates, and one has to deal with the full complexity of 3+1 numerical 
relativity. 

The fundamental ideas of the numerical solution of the conformal field equa- 
tions have been laid out by Frauendiener in this volume, and in a Living Review 
||, and he has also discussed his code to treat spacetimes with a hypersurface- 
orthogonal Killing vector and toroidal J^'s. The purpose of the present article 
is to show the status of numerical simulations based on the conformal field 
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equations in 3D - i.e. three space dimensions without assuming any continuous 
symmetries - and to discuss what is needed in order to render this approach a 
practical tool to investigate physically interesting spacetimes. By making future 
null infinity accessible to (completely regular and well defined) local computa- 
tions, the approach excels at the extraction of radiation - e.g. the quantities to 
hopefully be measured within the next years by new large-scale detectors. One of 
the main pedagogical goals will be to explain the challenges of numerical relativ- 
ity and to highlight some open problems related to constructing hyperboloidal 
initial data and actually carrying out long-time stable and accurate simulations. 
For an even more condensed account of the conformal approach to numerical 
relativity see [0. 

The organization is as follows: Sec. || introduces the algorithms developed 
in the last years by Peter Hiibner (the radiation extraction procedure, which 
I will only mention briefly, is based on work of Hiibner and Marsha Weaver), 
and implemented in a set of codes by Peter Hiibner (who has recently left the 
field). All results presented here have been obtained with these codes, which 
Hiibner has described in a series of articles H|],[l(|[Ll|. Sec. || will start with a 
brief description of the evolution of weak initial data which possess a regular 
point i + representing future timelike infinity, based on the work of Hiibner . 
Then I will discuss the evolution of slightly stronger initial data which exhibit 
various problems that will have to be solved, e.g. the choice of gauge, and use 
this as a starting point for discussing the main current problems. In Sec. [| purely 
computational aspects of this project will be discussed, and in Sec. ^ I will sum 
up the current status and sketch a possible roadmap for further work. 



2 Algorithms 

2.1 Problem Overview 

The conformal field equations are formulated in terms of an unphysical Lorentzian 
metric g ab defined on an unphysical manifold M. which gives rise to a physical 
metric g ab — ^~ 2 9ab, where the conformal factor £2 is determined by the equa- 
tions. The physical manifold M. is then given by A4 = {p € A4 \ f?(p) > 0}. 

Contrary to the formalism used by Fraucndiener in his contribution, we use 
a metric based formulation of the conformal field equations: 

V a i? bc - X7 b R ac = ((V a i?) g bc - (V b R) g ac ) ~ 2 (V d f2) d abc d , (1) 
y d d ab c d = 0, (2) 

V Q V b /2 = ^g ab V C V C Q - - R ab Q, (3) 

iv a (V 6 V b ^) =-^R ab V h Q~^QV a R-^V a f2R, (4) 

R a bc d = f?d abc d + (g ca Rb d — g C bR a d — g d a Rbc + g d bR a c)/2 

R 

+(gc a gb d - g C bg a d )j^, (5) 



Problems and Successes in the Numerical Approach 



3 



6J?V Q V a l2 = 12(V Q I2) (V a ft) - ft 2 R. (6) 

Here the Ricci scalar R of g a b is considered a given function of the coordinates. 
For any solution (g a b, Rab, d a bc d , ft), Rab is the traceless part of the Ricci tensor, 
and fl d a bc d the Weyl tensor of g a b- Note that the equations are regular even for 
ft = 0. These "conformal field equations" render possible studies of the global 
structure of spacetimes, e.g. reading off radiation at null infinity, by solving 
regular equations. 

The 3+1 decomposition of the conformal geometry can be carried out as 
usual in general relativity, e.g. 

9ab = h ab - n a n b = ft 2 (h ab - h a h b ), 

where h a b and h a b are the Riemannian 3-metrics induced by g a b respectively 
g a b on a spacelike hypersurface with unit normals n a , and equivalently h a = 
ftn a (our signature is (—,+,+,+)). The relation of the extrinsic curvatures 
(k ab = \Cnhab, kab = \C-nKb) is then easily derived as k ab = ft(k ab + ftoKb), 
where fto = n a V a /2. Note that for regular components of h a b and k a b, the 
corresponding components of h a b and k a b with respect to the same coordinate 
system will in general diverge due to the compactification effect. However for 
the coordinate independent traces k = h ab k a b, k — h ab k a b we get 

Qk = (k + 3ft ), 

which can be assumed regular everywhere. Note that at ^f, k = — 3fto- Since 
is an ingoing null surface (with (V a i7)(V a ]7) = but V a J? ^ ), we have 
that fto < at J? + . It follows that k > at . We will thus call regular 
spacelike hypersurfaces in M hyperboloidal hypersurfaces, since in M. they are 
analogous to the standard hyperboloids t 2 ~ x 2 ~ y 2 — z 2 = k 2 in Minkowski space, 
which provide the standard example. Since such hypersurfaces cross J? but are 
everywhere spacelike in A^, they allow to access and radiation quantities 
defined there by solving a Cauchy problem (in contrast to a characteristic initial 
value problem which utilizes a null surface slicing). Note that hyperboloidal 
hypersurfaces which cross ^ + are only Cauchy surfaces for the future domain 
of dependence of the initial slice of M, we therefore call our studies semiglobal. 

We will not discuss the full 3 + 1 equations here for brevity, but rather 
refer to What is important, is that the equations split into symmetric hyper- 
bolic evolution equations plus constraints which are propagated by the evolution 
equations The evolution variables are h a b, k a b, the connection coefficients 
7 a 6c projections (01) ^ a = n b h a c Rb c and (1,1) R a b = h a c h d Rbd of the traceless 
4-dimensional Ricci tensor Rbd, the electric and magnetic components of the 
rescaled Weyl tensor d a bc d E a b, B a b, as well as ft, fto, V a ft, V a V Q i? - in total 
this makes 57 quantities. In addition the gauge source functions q, R and N a 
have to be specified, in order to guarantee symmetric hyperbolicity they are given 
as functions of the coordinates. Here q determines the lapse as N = e^Vdet h 
and N a is the shift vector. The Ricci scalar R can be thought of as implicitly 
steering the conformal factor fl. 
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The numerical treatment of the constraints and evolution equations will be 
described below. But before, let us spend some time on general considerations 
about the treatment of null infinity. Since fl is an evolution variable and not 
specified a priori, J will in general not be aligned with grid points, and interpo- 
lation has to be used to evaluate computed quantities at locations of vanishing 
fl. For the physically interesting case of modeling an isolated system, "physical" 
.y - i.e. the component of that idealizes us outside observers and our gravita- 
tional wave detectors (neglecting cosmological effects such as redshift etc.) - has 
spherical topology. There may be more than one component of ^ ', i.e. ad diti onal 
spherical components associated with "topological black holes" (see Sec. 2.3). In 
principle it is possible of course to control the movement of J? through the grid 
by the gauge choice - see for how to achieve such J? fixing within Frauen- 
diencr's formulation. An example would be the so called freezing, where 
does not change its coordinate location. 

What is the significance of how J? moves through the grid? This question is 
directly related to the question for the global structure of spacetime. Although 
many questions are left open, the present understanding of the global structure of 
generic vacuum spacetimes, which can be constructed from regular initial data, 
does provide some hints. First, note that spacetimes which are asymptotically 
flat in spacelike and null directions - i.e. isolated systems - do not necessarily 
have to be asymptotically flat in timelike directions. An example would be a 
spacetime that contains a star or a black hole. In such cases where the end 
state of the system is not flat space, we can not expect the conformal spacetime 
(M.,g a b) to contain a regular point i + . In the case of sufficiently weak data 
however, Friedrich has shown in |L3| that a regular point i + will exist - consistent 
with our intuition. The global structure is then similar to Minkowski space. 
The standard conformal compactification of Minkowski space is discussed in 
textbooks (see e.g. ||) as a mapping to the Einstein static universe. There J" 
moves inward and contracts to a point within finite coordinate time. In order to 
resolve such situations it seems most appropriate to choose a gauge which mimics 
this behavior, i.e. where J? contracts to a point after finite coordinate time. The 
boundary of the computational domain is set in the unphysical region and the 
physical region contracts in coordinate space. Accordingly, the initial data are 
also extended beyond the physical region of spacetime. It is this scenario which 
is best understood so far, and which is presented in some more detail in Sec. |3j. 

For sufficiently strong regular data it is known that singularities develop fll4f 
- according to the cosmic censorship conjecture such singularities should 
generically be hidden inside of black holes. For such data we cannot expect a 
regular point i + to exist. In the case when i + is singular (and not much else is 
currently known even about the i + of Kruskal spacetime - see however Bernd 
Schmidt's contribution in this volume) we have to expect structure like sharp 
gradients near i + , which makes it unlikely that we can afford to significantly 
reduce the size of the physical region in coordinate space (at least not without 
adaptive mesh refinement - a technology not yet available for 3D evolutions). 
A J^-freezing gauge may be appropriate for such a situation. Furthermore, phe- 
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nomcna like quasi-normal ringdown, or the orbital motion of a two-black hole 
system suggest that the numerical time-coordinate better be adapted to the 
intrinsic time scale of the system. Associated with quasi-normal ringdown for 
example is a fixed period in Bondi-time, which suggests Bondi time as a time 
coordinate near / in a situation dominated by ringdown. Thus, for black hole 
spacetimes it might turn out that the best choice of gauge fixes J? to a particular 
coordinate position, and shifts i + into the infinite future []. It could be possible 
that in such a case the boundary can be chosen to either coincide with J^, or to 
be put just a small number of gridpoints outside, which would raise the ques- 
tion for an evolution algorithm that does not require a topologically rectangular 
grid. We see that the optimal choice of numerical algorithms and gauges may be 
tightly related to the global structure of the investigated spacetimes - which is 
actually one of the main questions our simulations should be able to answer! 

2.2 Construction of "extended" hyperboloidal initial data 

Evolution of a solution to the Einstein equations starts with a solution to the 
constraints. The constraints of the conformal field equations (see Eq. (14) of Ref. 
[§]) are regular equations on the whole conformal spacetime (Ai,g a b)- However, 
they have not yet been cast into a standard type of PDE system, such as a 
system of elliptic PDEs One therefore resorts to a 3-step method jlO| : 

1. Obtain data for the Einstein equations: the first and second fundamental 
forms h a b and k a b induced on S by g a b, corresponding in the compactified 
picture to h a b, k a b and fl and Hq. This yields so-called "minimal data". 

2. Complete the minimal data on £ to data for all variables using the conformal 
constraints - in principle this is mere algebra and differentiation. 

3. Extend the data from £ to £ in some ad hoc but sufficiently smooth and 
"well-behaved" way. 

In order to simplify the first step, the implementation of the code is restricted 
to a subclass of hyperboloidal slices where initially k a b is pure trace, fc a {, = ^h a bk. 
The momentum constraint V h k a b — V a /c = then implies k — const. ^ 0. We 
always set k > 0. In order to reduce the Hamiltonian constraint 

(3) fl + k 2 = k ab k ab 

to one elliptic equation of second order, we use a modified Lichnerowicz ansatz 

h a b = Q~ 2 4> h a b 

with two conformal factors fl and <p. The principal idea is to choose h a b and i?, 
and solve for 0, as we will describe now. First, the "boundary defining" function 
I? is chosen to vanish on a 2-surface S - the boundary of £ and initial cut of J? 

1 Note however, that in the absence of results, we are left to speculation here. 

2 work toward this goal is reported in this volume by Adrian Butscher. 
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- with non- vanishing gradient on S. The topology of S is chosen as spherical for 
asymptotically Minkowski spacetimes. Then we choose h a b to be a Ricmannian 
metric on S, with the only restriction that the extrinsic 2-curvature induced by 
h a b on S is pure trace, which is required as a smoothness condition With 
this ansatz h a b is singular at S, indicating that S represents an infinity. The 
Hamiltonian constraint then reduces to the Yamabe equation for the conformal 
factor <f>: 

4 /2 2 V a V Q 0-4 i7(V a J7)(V a 0)- Q (3) R Q 2 + 2QAQ - 3(V a J?)(V a r2)^ = \k 2 4> b . 

This is a semilinear elliptic equation - except at S, where the principal part 
vanishes for a regular solution. This however determines the boundary values as 

4 = 9fc- 2 (V a J7)(V a J7). (7) 

Existence and uniqueness of a positive solution to the Yamabe equation and the 
corresponding existence and uniqueness of regular data for the conformal field 
equations using the approach outlined above have been proven by Andersson, 
Chrusciel and Friedrich fl6| . Solutions to the Yamabe equation - and thus min- 
imal initial data - can either be taken from exact solutions or from numerical 
solutions of the Yamabe equation. Exact solutions which possess a J? of spherical 
topology have been implemented for Minkowski space and for Kruskal spacetime 

- see the contribution of Bernd Schmidt in this volume. These solutions are de- 
fined even outside of and thus can directly be completed to initial data for 
all variables by using the conformal constraints. 

If the Yamabe equation is solved numerically, the boundary has to be chosen 
at S, the initial cut of with boundary values satisfying Eq. ([?]). If the equation 
would be solved on a larger (more convenient Cartesian) grid, generic boundary 
conditions would cause the solution to lack sufficient differentiability at S, see 
Hiibner's discussion in Wfa. This problem is due to the degeneracy of the Yamabe 
equation at S. Unfortunately, this means that we have to solve an elliptic problem 
with spherical boundary. This problem is solved by combining the use of spher- 
ical coordinates with pseudo-spectral collocation methods. In pseudo-spectral 
methods the solution is expanded in (analytically known) basis functions - here 
a Fourier series for the angles and a Chebychev series for the radial coordinate. 
For an introduction to pseudo-spectral methods see e.g. [|l7[jl8| , |l^] . This allows 
to take care of coordinate singularities in a clean way provided that all ten- 
sor components are computed with respect to a regular (e.g. Cartesian) basis 
and that no collocation points align with the coordinate singularities. Another 
significant advantage of spectral methods is their fast convergence: for smooth 
solutions they typically converge exponentially with resolution. The necessary 
conversions between the collocation and spectral representations are carried out 
as fast Fourier transformations with the FFTW library [^0| . The nonlinearities 
are dealt with by a Newton iteration pl| , the resulting linear equations are 
solved by an algebraic multigrid linear solver (the AMG library |p2fl ). 

The constraints needed to complete minimal initial data to data for all evolu- 
tion variables split into two groups: those that require divisions by the conformal 
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factor fl to solve for the unknown variable, and those which do not. The latter 
do not cause any problems and can be solved without taking special care at 
Q = 0. The first group, needed to compute {1,1) R, E a b and B a b, however does 
require special numerical techniques to carry out the division, and furthermore 
it is not known if their solution outside of ^ actually allows solutions which are 
sufficiently smooth beyond J 1 ' , Thus, at least for these we have to find some ad- 
hoc extension. Note that in the case of analytical minimal data, the additional 
constraints are solved on the whole time evolution grid. 

The simplest approach to the division by fl would be an implementation of 
l'Hospital's rule, however this leads to nonsmooth errors and consequently to a 
loss of convergence || . Instead Hiibner |J has developed a technique to replace 
a division g = f / fl by solving an elliptic equation of the type (actually some 
additional terms added for technical reasons are omitted here for simplicity) 

V a V a (fl 2 g - flf) = 

for g. For the boundary values f2 2 g — flf = 0, the unique solution is g = j j fl. 
For technical details see |LC|| . The resulting linear elliptic equations for g are 
solved by the same numerical techniques as the Yamabe equation. For technical 
details see Hiibner Q. 

Finally, we have to extend the initial data to the full Cartesian spatial grid 
in some way. Since solving all constraints also outside of J 1 will in general not be 
possible in a sufficiently smooth way, we have to find an ad hoc extension, which 
violates the constraints outside of J? but is sufficiently well behaved to serve as 
initial data. The resulting constraint violation is not necessarily harmful for the 
evolution, since J? causally disconnects the physical region from the region of 
constraint violation. On the numerical level, errors from the constraint violating 
region will in general propagate into the physical region, but if our scheme is 
consistent, these errors have to converge to zero with the convergence order of 
the numerical scheme (fourth order in our case). There may still be practical 
problems, that prevent us from reaching this aim, of course: making the ad-hoc 
extension well behaved is actually quite difficult, the initial constraint violation 
may trigger constraint violating modes in the equations, which take us away 
from the true solution, singularities may form in the unphysical region, etc. 

Since the time evolution grid is Cartesian, its grid points will in general 
not coincide with the collocation points of the pseudo-spectral grid. Thus fast 
Fourier transformations can not be used for transformation to the time evolution 
grid. The current implementation instead uses standard discrete ("slow") Fourier 
transformations, which typically take up the major part of the computational 
effort of producing initial data. 

It turns out, that the combined procedure works reasonably well for certain 
data sets. For other data sets the division by fl is not yet solved in a satisfac- 
tory way, and constraint violations are of order unity for the highest available 
resolutions. In particular this concerns the constraint VbE a h — — {3) e a bck bd Bd c 
(Eq. (14d) in ||), since E a b is computed last in the hierarchy of variables and 
requires two divisions by fl. Further research is required to analyze the problems 
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and either improve the current implementation or apply alternative algorithms. 
Ultimately, it seems desirable to change the algorithm of obtaining initial data 
to a method that solves the conformal constraints directly and therefore does 
not suffer from the current problems. This approach may of course introduce 
new problems like an elliptic system too large to be handled in practice. 

2.3 Black hole initial data 

Since the standard definition of a black hole as the interior of an event horizon 
is a global concept, it is a priori not clear what one should consider as "black 
hole initial data" . In practice, the singularity theorems H] and the assumption 
of cosmic censorship ]l5| usually lead to the identification of "black hole initial 
data" with data that contain apparent horizons, and to associate the number of 
apparent horizons with the number of black holes in the initial data ^. 

A common strategy to produce apparent horizons is to use topologically 
nontrivial data, that is data which possess more than one asymptotically flat re- 
gion. In the time-symmetric case such data obviously possess a minimal surface! 
Asymptotic ends that extend to spatial infinity i° are relatively easy to produce 
by compactification methods, see e.g. ]27| , p8| . p9pc| ] or the contribution of Dain in 
this volume. From the numerical point of view it is important that the topology 
of the computational grid is independent of the number of asymptotic regions or 
apparent horizons considered: suitable regularization procedures allow to treat 
spatial infinities as grid points. 

In the current approach to the hyperboloidal initial value problem, where 
first the Yamabe equation needs to be solved, the grid topology does depend on 
the number of topological black holes - in this case the number of initial cuts of 
J^'s, which have spherical topology. One option would be of course to combine 
both ingredients and consider "mixed asymptotics" initial data, which extend 
to the physical .y and to unphysical interior spacelike infinities which only serve 
the purpose of acting as "topological sources" for apparent horizons. 

Another option, suggested by Hiibner in p0| , is to generalize the current 
code for the initial data, which only allows for one cut of y, which has spherical 
topology, to multiple ,fs of spherical topology. For the case of one black hole this 
would correspond to the relatively simple modification to S 2 x R topology. For the 
case of two black holes one could implement the Schwarz alternating procedure 
(as described in Sec. 6.4.1 of Ref. |24j]) to treat three ^'s with three coordinate 
patches, where each patch is adapted to a spherical coordinate system with its 
.y . A more practical approach (at least to get started) could be to produce 
topologically trivial black hole initial data. Since we expect physical black holes 
to result from the collapse of topologically trivial regular initial data, such data 
would in some sense be more physical. Theorems on the existence of apparent 
horizons in Cauchy data have been presented by Beig and O Murchadha in p7j . 
Numerical studies in this spirit have been performed by the author |2{|. Such 

3 Note that both the number of components of the event and apparent horizon are 
slice-dependent . 
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data could in principle be produced with the current code once it gets coupled 
to an apparent horizon finder. For the hyperboloidal initial value problem it is 
actually not known, whether such data actually exist, but it seems physically 
reasonable. Finding such data numerically by parameter studies would be an 
interesting result in itself. A natural question in this context is whether there 
is any qualitative difference between "topological" and "non-topological" black 
holes outside of the event horizon, e.g. regarding their waveforms? 

2.4 Numerical setup for evolutions 

The time evolution algorithm is an implementation of a standard fourth order 
method of lines (see e.g. (2j|), with centered spatial differences and Runge- 
Kutta time integration. Additionally, a dissipation term of the type discussed 
in theorems 6.7.1 and 6.7.2 of Gustafsson, Kreiss and Oliger J|5| is added to 
the right-hand-sides to damp out high frequency oscillations and keep the code 
numerically stable. Numerical experiments show that usually small amounts if 
dissipation are sufficient (the dissipation term used contains a free parameter), 
and do not change the results in any significant manner. 

A particularly subtle part of the evolution usually is the boundary treatment. 
In the conformal approach we are in the situation that the boundary is actually 
situated outside of the physical region of the grid - this is one of its essential 
advantages! In typical explicit time evolution algorithms, such as our Runge- 
Kutta method of lines, the numerical propagation speed is actually larger than 
the speed of all the characteristics (in our case the speed of light). Thus does 
not shield the physical region from the influence of the boundary - but this 
influence has to converge to zero with the convergence order of the algorithm 
- fourth order in our case. One therefore does not have to choose a "physical" 
boundary condition, the only requirements are stability and "practicality" - 
e.g. the boundary condition should avoid, if possible, the development of large 
gradients in the unphysical region to reduce the numerical "spill over" into the 
physical region, or even code crashes. 

The current implementation relies on a "transition layer" in the unphysical 
region, which is used to transform the rescaled Einstein equations to trivial 
evolution equations, which are stable with a trivial copy operation at outermost 
gridpoint as a boundary condition (see Ref. 0] for details and references). We 
thus modify the evolution equations according to replace 

dtf + A%f -b = Q d t f + a{Q) (A%f - b) = 0, 

where a is chosen as a(J?) = for SI < SIq < Sl\ < and 1 for fl > fl\. This 
procedure works reasonably well for weak data, however there are some open 
problems. One is, that the region of large constraint violations outside of may 
trigger constraint violating modes of the equations that can grow exponentially. 
Another problem ist that a "thin" transition zone causes large gradients in the 
coefficients of the equations - thus eventually leading to large gradients in the 
solution, while a "thick" transition zone means to loose many gridpoints. If no 
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transition zone is used at all, and the Cartesian grid boundary touches the 
ratio of the number of grid points in the physical region versus the number of 
grid points in the physical region is already 7r/6 « 0.52. 

2.5 Physics Extraction 

Extracting physics from a numerical solution to the Einstein equations is a non- 
trivial task. Results typically show a combination of physics and coordinate 
effects which are hard to disentangle, in particular in the absence of a back- 
ground geometry or preferred coordinate system. In order to understand what 
is going on in a simulation, e.g. to find "hot spots" of inaccuracy or instability 
or bugs in an algorithm, it is often very important to visualize the "raw" data 
of a calculation. Here the visualization of scalar and in particular tensor fields in 
3D is a subtle task in itself. But beyond that, one also wants ways to factor out 
coordinate effects in some way, and ideally access physical information directly. 

One way commonly used to partially factor our coordinate effects is to look at 
curvature invariants, another possibility is to trace geodesies through spacetime. 
In the current code this is done by concurrent integration of geodesies by means 
of the same 4th order Runge-Kutta scheme used already in the method of lines. 
Both null and timelike geodesies, as well as geodesies of the physical and rescaled 
metrics can be computed, and various quantities such as curvature invariants are 
computed by interpolation along the geodesies. 

Particularly important are null geodesies propagating along J^, since they 
can be used to define a Bondi system and thus compute radiation quantities 
such as the Bondi mass or news. Note that the foliation of spacetime cho- 
sen for evolution will in general not reproduce cuts of of constant Bondi 
time. Hiibner has therefore implemented postprocessing algorithms (using the 
IDL programming language/software system) which construct slices of constant 
Bondi time in the data corresponding to the null geodesies propagating on 
by interpolation (the algorithms are based on unpublished work of Hiibner and 
Weaver). This evolution of geodesies is illustrated by Fig. |l|, which shows three 
timelike geodesies originating with different initial velocities at the same point 
(x , y , z ) = (^=, j^=, j^=) meeting a generator of ^ at i + . 

3 Results for weak data 

In this section I will discuss results of 3D calculations for initial data which 
evolve into a regular point i + , and which thus could be called "weak data". 
Bcrnd Schmidt presents results for the Kruskal spacetime in this volume (see 
also j26j). The initial conformal metric is chosen in Cartesian coordinates as 

ds 2 = f 1 + jO 2 (x 2 + 2y 2 )^j dx 2 + dy 2 + dz 2 . (8) 

The boundary defining function is chosen as J? = (l — (x 2 + y 2 + z 2 )) /2, it is 
used to satisfy the smoothness condition for the conformal metric at J* . 
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1.2 1.4 



Fig. 1. Three timelike geodesies ol different initial velocities (starting at x = ^jjg) 
meet a generator of J? (starting at x — 1) at future timelike infinity i + . 



These data have been evolved previously by Hiibner for A = 1 as reported 
in pjj . For the gauge source functions he has made the "trivial" choice: R = 0, 
N a = 0, q = 0, i.e. the conformal spacetime has vanishing scalar curvature, 
the shift vanishes and the lapse is given by N = e q Vdeth = ydet h. This 
simplest choice of gauge is completely sufficient for A = 1 data, and has lead 
to a milestone result of the conformal approach - the evolution of weak data 
which evolve into a regular point i + of M., which is resolved as a single grid 
cell. With this result Hiibner has illustrated a theorem by Friedrich, who has 
shown that for sufficiently weak initial data there exists a regular point i + of M. 
[[l3). The complete future of (the physical part of) the initial slice can thus be 
reconstructed in a finite number of computational time steps. This calculation is 
an example of a situation for which the usage of the conformal field equations is 
ideally suited: main difficulties of the problem are directly addressed and solved 
by using the conformal field equations. 

The natural next question to ask is: what happens if one increases the ampli- 
tude A? To answer this question, I have performed and analyzed runs for integer 
values of A up to A = 20. The results presented here have been produced with 
low resolutions of 80 3 (but for higher or slightly lower resolutions we essentially 
get the same results) . For convergence tests of the code see P,|l0| 
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While for 

A = 1,2 the code continues beyond i + without problems, for all higher ampli- 
tudes the "trivial" gauge leads to code crashes before reaching i + . Here by "code 
crash" we mean that computational values get undefined, e.g. the code produces 
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"not a number" (NaN) values. While the physical data still decay quickly in 
time, a sharp peak of the lapse develops outside of and crashes the code after 
Bondi time - 8(320M) for A = 3 and - 1.5 (3 A/) for A = 20 (here M is the ini- 
tial Bondi mass). In Figs. || and |3] the lapse N is plotted for runs with A = 1 and 
A = 5. While for A = 1 the lapse only shows significant growth after t = 1 (i + 
is located at x = y = z = 0, t = 1), for A = 5 a very sharp peak grows outside 
of J 2 " and crashes the code at t = 0.9076. Where does this rapid growth come 
from? Note that the initial conformal metric Eq. (|J) shows significant growth 
outside of J 1 . Combined with the lapse N = V det h this leads to a growth of the 
lapse toward the grid boundary. In the present case a positive feedback effect 
with a growth of metric components in time seems to be responsible for the 
eventual crash of the code. Note that this feedback only takes place in a small 
region outside of J 1 - further outward it is prevented by the transition to trivial 
evolution equations. 

Fig. ^ shows a cure of the problem: a modified gauge source function q = 
—r 2 1 a (N = e~ r l a \J 'det h) with a — 1 leads to a very smooth lapse (and cor- 
respondingly also to smooth metric components). Note that in Fig. [|, due to 
the different lapse, the point i + is not located at t = 1. The value of a = 1 
here is found by moderate tuning of a to a best value (significantly decreasing 
or increasing a crashes the code before i + is reached) . Unfortunately, this modi- 
fication of the lapse is not sufficient to achieve much higher amplitudes. As A is 
increased, the parameter a requires more fine tuning, which was only achieved 
for A < 8. For higher amplitudes the code crashes with significant differences in 
the maximal and minimal Bondi time achieved, while the radiation still decays 
very rapidly and the news scales almost linearly. Furthermore, the curvature 
quantities do not show excessive growth - it is thus natural to assume that we 
are still in the weak-field regime, and the crash is not connected to the formation 
of an apparent horizon or singularity. 

These results suggest that in order to model a gauge source function q that 
would allow to evolve up to i + , one would need more that one parameter, e.g. at 
least 3 parameters for a non-isotropic ansatz such as q — — (x 2 / a 2 +y 2 / b 2 + z 2 / c 2 ) 
or something similar. To tune 3 or more parameters for each evolution seems how- 
ever computationally prohibitive. While some improvement is obviously possible 
through simple non-trivial models for the lapse (or other gauge source functions), 
this approach seems very limited and more understanding will be necessary to 
find practicable gauges. An interesting line of research would be to follow the 
lines of Ref. |52| in order to find evolution equations for the gauge source func- 
tions which avoid the development of pathologies. A particular aim would be to 
find equations such that the resulting system of evolution equations is symmetric 
hyperbolic. 

Fig. H shows the news function from three different runs: for A = 1,5 and 
q = 0, and for A = 5, q = —r 2 /a. The news from the runs for A = 1 have been 
multiplied by a factor of 25, which would exactly compensate scaling with the 
amplitude in the linear regime. We see that the three curves line up very well 
initially. The line for A = 5 and q = deviates significantly to larger values of 
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the news when the runs starts to get inaccurate, but at this time most of the 
physical radiation has already left the system. The curves from A = 1/q = 
and A = 5/q = —r 2 /a line up perfectly until the value of the news drops below 
10~ 6 , where the curves level off at different values, due to numerical inaccuracy. 

Fig. |^ shows the Bondi mass for this situation, again the A = 1 curve is scaled 
by a factor of 25: Again we see the quick decay of a sharp pulse of radiation. 
There is no particular structure except falloff at late times, the deviation of the 
curves at late times seems to be caused by numerical inaccuracy, in particular 
in the computation of the Bondi mass. 



Fig. 2. Lapse N for the x-axis, for amplitude A — 1 and R = N a = q = 0. The 
contour line marks {Q — 0). 



4 Computational aspects 

In this section I will give a brief description of some computational aspects, 
such as the computational resources needed to carry out simulations in 3 spatial 
dimensions. Computations of this scale rely on parallel processing, that means 
execution of our algorithms is spread over different CPU's. From a simplistic 
point of view there are two ways to program for parallel execution: we only take 
care of parallelizing the algorithm - but we require that all CPU's can access the 
same memory - or we both parallelize the algorithm and the data structures, 
and separate the total data into smaller chunks that fit into the local memory of 
each processor. The first alternative requires so called shared memory machines, 
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Fig. 3. Lapse N for the :r-axis, for A = 5 and R = N a = q = 0. The contour line 
marks J> {Q = 0). 




Fig. 4. Lapse AT for the a;-axis, for A — 5 and R — N a = 0, a = — r 2 . The contour 
line marks {Q = 0). 
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where the operating system and hardware take care of making data accessible 
to the CPU's consistently, taking care of several layers of main and cache mem- 
ory (which gets increasingly difficult and expensive as the size of the machine 
is increased). The present code has been implemented using a shared memory 
programming model. The advantage is that this can generally be somewhat eas- 
ier to program, and avoids overheads in memory. The disadvantage is the high 
cost of such systems, which makes them difficult to afford and thus nonstandard 
for most large academic parallel applications. The second alternative, usually 
referred to as distributed memory, requires more work to be done by the pro- 
grammer, but more flexibly adapts to different kinds of machines such as clusters 
of cheap workstations commonly available in academic environments. While this 
approach usually implies a larger overhead in total memory requirements, speed 
and programming complexity, it is currently the only approach capable of scaling 
from small to very large simulations. For a general introduction to the issues of 
high performance computing, Ref. |57j provides a good starting point. 

So, how much memory do we need? Let us assume a run with 150 3 grid points 
(the size of the largest simulations carried out with the present code so far). The 
current implementation of the fourth order Runge-Kutta algorithm uses 4 time 
levels and a minimum of 62 grid functions (57 variables and 5 gauge source func- 
tions). In double precision this amounts to ~ 6.7 GByte. Temporary variables, 
information on geodesies and various overheads result in a typical increase of 
memory requirements by roughly 20%. For 150 time steps (approximately what 
it takes to reach i + for weak data) the total amount of processed data then 
corresponds to roughly 1 Terabyte! 

If we half the grid spacing, the allocated memory increases by a factor of 
8 (neglecting overheads), the total amount of processed data by a factor of 
16, and the total required CPU time also by a factor of 16, while the error 
reduces by a factor of 16 - if we are already in the convergent regime! Given 
that the biggest academic shared memory machines in Germany have 16 GByte 
of memory available (the AEI's Origin 2000 and the Hitachi SR8000 at LRZ 
in Munich) this shows that the margin for increasing resolution is currently 
quite small. Such an increase in resolution will however be necessary to resolve 
physically interesting situations with more structure, such as a black hole, or 
two merging black holes. A move toward distributed memory processing will 
therefore be likely in the long run. 

The current software-standard for distributed (scalable) computing is MPI 
(message passing interface) Q . Unfortunately, writing large scale sophisticated 
codes in MPI is very time consuming. However, several software packages are 
available which introduce a software layer between the application programmer 
and MPI, and thus significantly reduce the effort to write parallel applications. 
Two prime examples are the Cactus Computational Toolkit |38| and PETSc 
[ B9) . While PETSc is a general purpose tool developed at Argonne National 
Laboratory as parallel framework for numerical computations, Cactus has been 
developed at the Albert Einstein Institute with numerical relativity in mind. 
While PETSc offers more support for numerical algorithms, in particular for 
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parallel elliptic solvers, Cactus already contains some general numerical relativity 
functionality like apparent horizon finders - but no support for generic numerical 
algorithms. Apart from its numerical relativity flavor, the Cactus computational 
toolkit also has the advantage of broad support for parallel I/O and large scale 
3D visualization. The ability to successfully mine tens or hundreds of gigabytes 
of data for relevant features is paramount to successful simulations in 3D. 

Among the essential problems in writing and maintaining large scientific 
codes are the software engineering aspects and the control of complexity. In other 
words, codes should be reasonably documented and maintainable. For large sci- 
entific codes written and maintained by part-time-programmer scientists this 
poses a significant challenge. Writing a clear, modular code that can be under- 
stood, maintained and extended to suit new scientific needs requires a good deal 
of design and planning ahead. For an introduction to software engineering issues 
see e.g. ^l[, ^2|. Another important issue for scientific codes is flexibility. Being 
able to do good science often depends on the ability to easily change algorithms, 
equations, discretization schemes etc. - without having to restructure the code, 
without high risk of introducing new bugs. In the present case examples for the 
need of trying different things would be experiments with different evolution 
equations (e.g. metric versus frame formalism), different boundary treatments 
or different elliptic solvers. 

5 Discussion 

The 3D numerical simulations performed so far show that the evolutions are 
numerically stable and quite robust. However, one of the main problems in nu- 
merical relativity is the stability of the constraint propagation: while the con- 
straints do propagate when they are satisfied identically initially, this assumption 
does not hold for numerical simulations. On the contrary - it seems to be quite 
typical observation that the constraints diverge exponentially, if the evolution 
does not start at the constraint surface. Preliminary results exhibit this behav- 
ior of resolution-independent exponential growth associated with a violation of 
the constraints also for the conformal approach. One of the major goals for the 
future thus has to be the improvement of the understanding of the constraint 
propagation equations, and an according modification of the evolution equations 
(see Ref. |33| for previous work in this direction). This is essentially an analytical 
problem, but will certainly require the numerical testing of ideas. 

Another area where new developments are necessary on the analytical side 
- along with numerical testing - is the problem of finding gauges that prevent 
pathologies like unnecessarily strong gradients. Ideally one would want to keep 
the symmetric hyperbolic character of the evolution system while allowing for a 
maximum of flexibility in writing down evolution equations for the gauge source 
functions. 

Well-posedness of the evolution equations is important - but by far not suf- 
ficient for numerical purposes. While well-posedness unfortunately still has not 
yet been shown rigorously for many formulations used in numerical relativity, 
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another important task seems to be to improve the understanding of the non- 
principal part of the equations, including their nonlinearities, in order to be able 
to construct numerically well-conditioned algorithms. 

The third area where significant progress seems necessary on the analytical 
side is the construction of initial data. Problems with the current algorithm 
which necessitates divisions by zero and an ad-hoc extension beyond have 
not yet been resolved. A possible road toward resolving these problems has been 
outlined by Butscher in this volume. 

An important role in improving the analytical understanding and in setting 
up numerical experiments will be played by the utilization of simplifications. 
Particularly important are spacetime symmetries and perturbative studies. A 
particularly interesting case to be studied actually is Minkowski space. Besides 
being an important case for code testing, it is used in current investigations to 
learn more about gauges and the stability of constraint propagation. 

More complicated are general spherically symmetric spacetimes. In the vac- 
uum case, this only leaves the Kruskal spacetime aside from Minkowski space 
- but understanding the gauge problem for Kruskal spacetime is an important 
milestone toward long-time black hole simulations. Moreover, spherical symme- 
try provides a natural testing ground for all kinds of new ideas, e.g. of how 
to treat the appearance of singularities, of how to treat the unphysical region, 
numerical methods, etc. 

An alternative route to simplification, which has been very successful in nu- 
merical relativity, is perturbative analysis, e.g. with Minkowski or Schwarzschild 
backgrounds. In the context of compactification this has been carried out nu- 
merically with characteristic codes in p5p^ ] , some of the problems that showed 
up there are likely to be relevant also for the conformal approach. 

What can we expect from the conformal approach in terms of physics results? 
Where can we expect contributions to our understanding of general relativity? 
One of the most important features of the conformal approach is that it excels 
at radiation extraction without ambiguities, and - at least in principle - enables 
numerical codes to study the global structure of spacetimes describing isolated 
systems. As has been demonstrated in this paper, in some weak field regime 
the code works well with relatively simple choices of gauge, and could be used 
to investigate some of the above problems. It could also provide a very clean 
way to study nonlinear deviations from linear predictions. For strong fields, in 
particular one or two black holes, the problem is much more difficult. An even 
more difficult problem is the investigation of the structure of singularities. In 
the spherically symmetric case this has been achieved by Hiibner |34| , but it is 
not clear whether these methods can be carried over to the generic case without 
symmetries, where the structure of the singularity has to be expected to be much 
more complicated. 

What is the roadmap for the future? As far as 3D simulations are concerned, 
I believe that one should try to go from relatively well controlled weak data 
to stronger data and try to identify and solve problems as they come up. In 
parallel, it will be important to study simplified situations, like spacetimes with 
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symmetries or linear perturbations with a mixture of analytical and numeri- 
cal techniques. Both lines of research are hoped to improve our understanding 
of issues associated with choosing the gauge source functions, and controlling 
the growth of constraints. Future 3D codes, aimed at producing novel physical 
results will also require a significant effort devoted to "computational engineer- 
ing", since flexible and solidly written codes are an absolute necessity for good 
computational science! These well known problems plaguing 3D numerical rela- 
tivity will have to be addressed and solved in the conformal approach in order to 
harvest its benefits. Developing the conformal approach to numerical relativity 
into a mature tool poses an important challenge for mathematical relativity: not 
only is the problem hard and requires long-term investments, it also requires 
to merge sophisticated mathematical analysis with computational engineering. 
The aim is to produce a solid handle on exciting new physics - and some of the 
physics will even be accessible to experiments. 
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